﻿#region << 版 本 注 释 >>
/*----------------------------------------------------------------
 * 版权所有 (c) 2022 北京超维景生物科技有限公司 保留所有权利。
 * 
 * 创建者：huangyang
 * 电子邮箱：huangyang@tvscope.cn
 * 创建时间：2023/2/16 10:42:19
 * 版本：V1.0.0
 * 描述：
 *
 * ----------------------------------------------------------------
 * 修改人：
 * 时间：
 * 修改说明：
 *
 * 版本：V1.0.1
 *----------------------------------------------------------------*/
#endregion << 版 本 注 释 >>

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;

namespace ImageK.Measure
{
    /** Curve fitting class based on the Simplex method in the Minimizer class
     *
     *
     *	Notes on fitting polynomial functions:
     *	(i) The range of x values should not be too far from 0, especially for higher-order polynomials.
     *	For polynomials of 4th order, the average x value should fulfill |xMean| < 2*(xMax-xMin).
     *	For polynomials of 5th order or higher, the x range should encompass x=0; and for
     *	7th and 8th order it is desirable to have x=0 near the center of the x range.
     *	(ii) In contrast to some fitting algorithms using the normal equations and matrix inversion, the
     *	simplex algorithm used here can cope with parameters having very different orders of magnitude,
     *	as long as the coefficients of the polynomial are within a reasonable range (say, 1e-80 to 1e80).
     *	Thus, it is usually not needed to scale the x values, even for high-order polynomials.
     *
     *	Version history:
     *
     *	2008-01-21: Modified to do Gaussian fitting by Stefan Woerz (s.woerz at dkfz.de).
     *	2012-01-30: Modified for external Minimizer class and UserFunction fit by Michael Schmid.
     *			  - Also fixed incorrect equation string for 'Gamma Variate' & 'Rodbard (NIH Image)',
     *			  - Added 'Inverse Rodbard', 'Exponential (linear regression)', 'Power (linear regression)'
     *				functions and polynomials of order 5-8. Added 'nicely' sorted list of types.
     *			  - Added absolute error for minimizer to avoid endless minimization if exact fit is possible.
     *			  - Added 'initialParamVariations' (order of magnitude if parameter variations) -
     *				this is important for safer and better convergence.
     *			  - Linear regression for up to 2 linear parameters, reduces the number of parameters handled
     *				by the simplex Minimizer and improves convergence.	These parameters can be an offset and
     *				either a linear slope or a factor that the full function is multiplied with.
     *	2012-10-07: added GAUSSIAN_NOOFFSET fit type
     *	2012-11-20: Bugfix: exception on Gaussian&Rodbard with initial params, bad initial params for Gaussian
     *  2013-09-24: Added "Exponential Recovery (no offset)" and "Chapman-Richards" (3-parameter;
     *              used e.g. to describe forest growth) fit types.
     *  2013-10-11: bugfixes, added setStatusAndEsc to show iterations and enable abort by ESC
     *  2015-03-26: bugfix, did not use linear regression for RODBARD
     *  2016-11-28: added static getNumParams methods
     *  2018-03-23: fixes NullPointerException for custom fit without initialParamVariations
     *  2018-07-19: added error function erf (=integral over Gaussian)
     *  2021-04-30: data points can have weights
     */
    public class CurveFitter:UserFunction
    {
        /** Constants for the built-in fit functions */
        public const int STRAIGHT_LINE = 0, POLY2 = 1, POLY3 = 2, POLY4 = 3,
            EXPONENTIAL = 4, POWER = 5, LOG = 6, RODBARD = 7, GAMMA_VARIATE = 8, LOG2 = 9,
            RODBARD2 = 10, EXP_WITH_OFFSET = 11, GAUSSIAN = 12, EXP_RECOVERY = 13,
            INV_RODBARD = 14, EXP_REGRESSION = 15, POWER_REGRESSION = 16,
            POLY5 = 17, POLY6 = 18, POLY7 = 19, POLY8 = 20,
            GAUSSIAN_NOOFFSET = 21, EXP_RECOVERY_NOOFFSET = 22, CHAPMAN = 23, ERF = 24;

        /** Nicer sequence of the built-in function types */
        public int[] sortedTypes = new int[]{ STRAIGHT_LINE, POLY2, POLY3, POLY4,
            POLY5, POLY6, POLY7, POLY8,
            POWER, POWER_REGRESSION,
            EXPONENTIAL, EXP_REGRESSION, EXP_WITH_OFFSET,
            EXP_RECOVERY, EXP_RECOVERY_NOOFFSET,
            LOG, LOG2,
            GAUSSIAN, GAUSSIAN_NOOFFSET, ERF,
            RODBARD, RODBARD2, INV_RODBARD,
            GAMMA_VARIATE,CHAPMAN
        };

        /** Names of the built-in fit functions */
        public static String[] fitList = {"Straight Line","2nd Degree Polynomial",
            "3rd Degree Polynomial", "4th Degree Polynomial","Exponential","Power",
            "Log","Rodbard", "Gamma Variate", "y = a+b*ln(x-c)","Rodbard (NIH Image)",
            "Exponential with Offset","Gaussian", "Exponential Recovery",
            "Inverse Rodbard", "Exponential (linear regression)", "Power (linear regression)",
            "5th Degree Polynomial","6th Degree Polynomial","7th Degree Polynomial","8th Degree Polynomial",
            "Gaussian (no offset)", "Exponential Recovery (no offset)",
            "Chapman-Richards", "Error Function"
        }; // fList, doFit(), getNumParams() and makeInitialParamsAndVariations() must also be updated

        /** Equations of the built-in fit functions */
        public static String[] fList = {
            "y = a+bx","y = a+bx+cx^2",									//STRAIGHT_LINE,POLY2
            "y = a+bx+cx^2+dx^3","y = a+bx+cx^2+dx^3+ex^4",
            "y = a*exp(bx)","y = a*x^b", "y = a*ln(bx)",				//EXPONENTIAL,POWER,LOG
            "y = d+(a-d)/(1+(x/c)^b)", "y = b*(x-a)^c*exp(-(x-a)/d)",	//RODBARD,GAMMA_VARIATE
            "y = a+b*ln(x-c)", "x = d+(a-d)/(1+(y/c)^b) [y = c*((x-a)/(d-x))^(1/b)]",  //LOG2,RODBARD2
            "y = a*exp(-bx) + c", "y = a + (b-a)*exp(-(x-c)*(x-c)/(2*d*d))", //EXP_WITH_OFFSET,GAUSSIAN
            "y = a*(1-exp(-b*x)) + c", "y = c*((x-a)/(d-x))^(1/b)",		//EXP_RECOVERY, INV_RODBARD
            "y = a*exp(bx)", "y = a*x^b",								//EXP_REGRESSION, POWER_REGRESSION
            "y = a+bx+cx^2+dx^3+ex^4+fx^5", "y = a+bx+cx^2+dx^3+ex^4+fx^5+gx^6",
            "y = a+bx+cx^2+dx^3+ex^4+fx^5+gx^6+hx^7", "y = a+bx+cx^2+dx^3+ex^4+fx^5+gx^6+hx^7+ix^8",
            "y = a*exp(-(x-b)*(x-b)/(2*c*c)))",						//GAUSSIAN_NOOFFSET
            "y = a*(1-exp(-b*x))",										//EXP_RECOVERY_NOOFFSET
            "y = a*(1-exp(-b*x))^c",									//CHAPMAN
            "y = a+b*erf((x-c)/d)"										//ERF; note that the c parameter is sqrt2 times the Gaussian
        };

        /** ImageJ Macro language code for the built-in functions */
        public static String[] fMacro = {
            "y = a+x*b","y = a+x*(b+x*c)",								//STRAIGHT_LINE,POLY2
            "y = a+x*(b+x*(c+x*d))","y = a+x*(b+x*(c+x*(d+x*e)))",
            "y = a*Math.exp(b*x)","y = a*Math.pow(x,b)", "y = a*Math.log(b*x)", //EXPONENTIAL,POWER,LOG
            "y = d+(a-d)/(1+Math.pow(x/c,b))", "y = b*Math.pow(x-a,c)*Math.exp(-(x-a)/d)",	//RODBARD,GAMMA_VARIATE
            "y = a+b*Math.log(x-c)", "y = c*Math.pow((x-a)/(d-x),1/b)",  	 //LOG2,RODBARD2
            "y = a*Math.exp(-b*x)+c", "y = a+(b-a)*Math.exp(-(x-c)*(x-c)/(2*d*d))", //EXP_WITH_OFFSET,GAUSSIAN
            "y = a*(1-Math.exp(-b*x))+c", "y = c*Math.pow((x-a)/(d-x),1/b)", //EXP_RECOVERY, INV_RODBARD
            "y = a*Math.exp(b*x)", "y = a*Math.pow(x,b)",						//EXP_REGRESSION, POWER_REGRESSION
            "y = a+x*(b+x*(c+x*(d+x*(e+x*f))))", "y = a+x*(b+x*(c+x*(d+x*(e+x*(f+x*g)))))",
            "y = a+x*(b+x*(c+x*(d+x*(e+x*(f+x*(g+x*h))))))", "y = a+x*(b+x*(c+x*(d+x*(e+x*(f+x*(g+x*(h+x*i)))))))",
            "y = a*Math.exp(-(x-b)*(x-b)/(2*c*c))",						//GAUSSIAN_NOOFFSET
            "y = a*(1-Math.exp(-b*x))",									//EXP_RECOVERY_NOOFFSET
            "y = a*Math.pow(1-Math.exp(-b*x),c)",						//CHAPMAN
            "y = a+b*Math.erf((x-c)/d)"									//ERF; note that the c parameter is sqrt2 times the Gaussian sigma
        };

		public double userFunction(double[] paras, double x)
        {
            throw new NotImplementedException();
        }

		/** Returns value of built-in 'fitType' formula value for parameters "p" at "x" */
		public static double f(int fitType, double[] p, double x)
		{
			switch (fitType)
			{
				case STRAIGHT_LINE:
					return p[0] + x * p[1];
				case POLY2:
					return p[0] + x * (p[1] + x * p[2]);
				case POLY3:
					return p[0] + x * (p[1] + x * (p[2] + x * p[3]));
				case POLY4:
					return p[0] + x * (p[1] + x * (p[2] + x * (p[3] + x * p[4])));
				case POLY5:
					return p[0] + x * (p[1] + x * (p[2] + x * (p[3] + x * (p[4] + x * p[5]))));
				case POLY6:
					return p[0] + x * (p[1] + x * (p[2] + x * (p[3] + x * (p[4] + x * (p[5] + x * p[6])))));
				case POLY7:
					return p[0] + x * (p[1] + x * (p[2] + x * (p[3] + x * (p[4] + x * (p[5] + x * (p[6] + x * p[7]))))));
				case POLY8:
					return p[0] + x * (p[1] + x * (p[2] + x * (p[3] + x * (p[4] + x * (p[5] + x * (p[6] + x * (p[7] + x * p[8])))))));

				//todo:
				// case EXPONENTIAL:
				// case EXP_REGRESSION:
				// 	return p[0] * Math.exp(p[1] * x);
				// case EXP_WITH_OFFSET:
				// 	return p[0] * Math.exp(-p[1] * x) + p[2];
				// case EXP_RECOVERY:
				// 	return p[0] * (1 - Math.exp(-p[1] * x)) + p[2];
				// case EXP_RECOVERY_NOOFFSET:
				// 	return p[0] * (1 - Math.exp(-p[1] * x));
				// case CHAPMAN:                               // a*(1-exp(-b*x))^c
				// 	double value = p[0] * (Math.pow((1 - (Math.exp(-p[1] * x))), p[2]));
				// 	//	Log.e("test", "values = " + value);
				// 	return value;
				// case GAUSSIAN:
				// 	return p[0] + (p[1] - p[0]) * Math.exp(-(x - p[2]) * (x - p[2]) / (2.0 * p[3] * p[3]));
				// case GAUSSIAN_INTERNAL:                     // replaces GAUSSIAN for the fitting process
				// 	return p[0] + p[1] * Math.exp(-(x - p[2]) * (x - p[2]) / (2.0 * p[3] * p[3]));
				// case GAUSSIAN_NOOFFSET:
				// 	return p[0] * Math.exp(-(x - p[1]) * (x - p[1]) / (2.0 * p[2] * p[2]));
				// case POWER:                                 // ax^b
				// case POWER_REGRESSION:
				// 	return p[0] * Math.pow(x, p[1]);
				// case LOG:
				// 	if (x == 0.0)
				// 		return -1000 * p[0];
				// 	return p[0] * Math.log(p[1] * x);
				// case RODBARD:
				// 	{                               // d+(a-d)/(1+(x/c)^b)
				// 		double ex = Math.pow(x / p[2], p[1]); //(x/c)^b
				// 		return p[3] + (p[0] - p[3]) / (1.0 + ex);
				// 	}
				// case RODBARD_INTERNAL:
				// 	{                   // d+a/(1+(x/c)^b) , replaces RODBARD of the fitting process
				// 		double ex = Math.pow(x / p[2], p[1]); //(x/c)^b
				// 		return p[3] + p[0] / (1.0 + ex);
				// 	}
				// case GAMMA_VARIATE:                         // b*(x-a)^c*exp(-(x-a)/d)
				// 	if (p[0] >= x) return 0.0;
				// 	if (p[1] <= 0) return Double.NaN;
				// 	if (p[2] <= 0) return Double.NaN;
				// 	if (p[3] <= 0) return Double.NaN;
				//
				// 	double pw = Math.pow((x - p[0]), p[2]);
				// 	double e = Math.exp((-(x - p[0])) / p[3]);
				// 	return p[1] * pw * e;
				// case LOG2:
				// 	double tmp = x - p[2];
				// 	if (tmp <= 0)
				// 		return Double.NaN;
				// 	return p[0] + p[1] * Math.log(tmp);
				// case INV_RODBARD:       // c*((x-a)/(d-x))^(1/b), the inverse Rodbard function
				// case RODBARD2:          // also used after the 'Rodbard NIH Image' fit
				// 	double y;
				// 	if (p[3] - x < 2 * Double.MIN_VALUE || x < p[0]) // avoid x>=d (singularity) and x<a (negative exponent)
				// 		y = fitType == INV_RODBARD ? Double.NaN : 0.0;
				// 	else
				// 	{
				// 		y = (x - p[0]) / (p[3] - x);        //(x-a)/(d-x) = ( (a-d)/(x-d) -1 )
				// 		y = Math.pow(y, 1.0 / p[1]);    //y=y**(1/b)
				// 		y = y * p[2];
				// 	}
				// 	return y;
				// case ERF:
				// 	return p[0] + p[1] * IJMath.erf((x - p[2]) / p[3]); //y=a+b*erf((x-c)/d)
				default:
					return 0.0;
			}
		}
	}
}
